﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Diagnostics;

namespace DeloneTriangulate
{
    namespace Triangula
    {
        class Vertex
        {
            public double x, y;

            protected Vertex() { }

            public Vertex(double x, double y)
            {
                this.x = x; this.y = y;
            }

            public double distance2To(Vertex other)
            {
                double dx = x - other.x;
                double dy = y - other.y;
                return dx * dx + dy * dy;
            }

            public double distanceTo(Vertex other)
            {
                return (double)Math.Sqrt(distance2To(other));
            }

            public override string ToString()
            {
                return string.Format("({0},{1})", x, y);
            }

            // бинарный плюс
            public static Vertex operator +(Vertex v1, Vertex v2)
            {
                return new Vertex(v1.x + v2.x, v1.y + v2.y);
            }

            // бинарный минус
            public static Vertex operator -(Vertex v1, Vertex v2)
            {
                return new Vertex(v1.x - v2.x, v1.y - v2.y);
            }

            // бинарное умножение
            public static Vertex operator *(double f, Vertex v)
            {
                return new Vertex(v.x * f, v.y * f);
            }

            // унарный плюс
            public static Vertex operator +(Vertex v)
            {
                return new Vertex(v.x, v.y);
            }

            // унарный минус
            public static Vertex operator -(Vertex v)
            {
                return new Vertex(-v.x, -v.y);
            }

            public static bool operator ==(Vertex v1, Vertex v2)
            {
                if ((v1.x == v2.x) && (v1.y == v2.y))
                    return true;
                return false;
            }

            public static bool operator !=(Vertex v1, Vertex v2)
            {
                if ((v1.x != v2.x) || (v1.y != v2.y))
                    return true;
                return false;
            }
        }

        class Triangulator
        {
            private List<Vertex> points;

            public Triangulator()
            {
            }

            private void Analyse(List<Vertex> suppliedPoints, Hull hull, List<Triad> triads, bool rejectDuplicatePoints, bool hullOnly)
            {
                if (suppliedPoints.Count < 3)
                    throw new ArgumentException("Number of points supplied must be >= 3");

                this.points = suppliedPoints;
                int nump = points.Count;

                double[] distance2ToCentre = new double[nump];
                int[] sortedIndices = new int[nump];

                // Choose first point as the seed
                for (int k = 0; k < nump; k++)
                {
                    distance2ToCentre[k] = points[0].distance2To(points[k]);
                    sortedIndices[k] = k;
                }

                // Sort by distance to seed point
                Array.Sort(distance2ToCentre, sortedIndices);

                // Duplicates are more efficiently rejected now we have sorted the vertices
                if (rejectDuplicatePoints)
                {
                    // Search backwards so each removal is independent of any other
                    for (int k = nump - 2; k >= 0; k--)
                    {
                        // If the points are identical then their distances will be the same,
                        // so they will be adjacent in the sorted list
                        if ((points[sortedIndices[k]].x == points[sortedIndices[k + 1]].x) &&
                            (points[sortedIndices[k]].y == points[sortedIndices[k + 1]].y))
                        {
                            // Duplicates are expected to be rare, so this is not particularly efficient
                            Array.Copy(sortedIndices, k + 2, sortedIndices, k + 1, nump - k - 2);
                            Array.Copy(distance2ToCentre, k + 2, distance2ToCentre, k + 1, nump - k - 2);
                            nump--;
                        }
                    }
                }

                Debug.WriteLine((points.Count - nump).ToString() + " duplicate points rejected");

                if (nump < 3)
                    throw new ArgumentException("Number of unique points supplied must be >= 3");

                int mid = -1;
                double romin2 = double.MaxValue, circumCentreX = 0, circumCentreY = 0;

                // Find the point which, with the first two points, creates the triangle with the smallest circumcircle
                Triad tri = new Triad(sortedIndices[0], sortedIndices[1], 2);
                for (int kc = 2; kc < nump; kc++)
                {
                    tri.c = sortedIndices[kc];
                    if (tri.FindCircumcirclePrecisely(points) && tri.circumcircleR2 < romin2)
                    {
                        mid = kc;
                        // Centre of the circumcentre of the seed triangle
                        romin2 = tri.circumcircleR2;
                        circumCentreX = tri.circumcircleX;
                        circumCentreY = tri.circumcircleY;
                    }
                    else if (romin2 * 4 < distance2ToCentre[kc])
                        break;
                }

                // Change the indices, if necessary, to make the 2th point produce the smallest circumcircle with the 0th and 1th
                if (mid != 2)
                {
                    int indexMid = sortedIndices[mid];
                    double distance2Mid = distance2ToCentre[mid];

                    Array.Copy(sortedIndices, 2, sortedIndices, 3, mid - 2);
                    Array.Copy(distance2ToCentre, 2, distance2ToCentre, 3, mid - 2);
                    sortedIndices[2] = indexMid;
                    distance2ToCentre[2] = distance2Mid;
                }

                // These three points are our seed triangle
                tri.c = sortedIndices[2];
                tri.MakeClockwise(points);
                tri.FindCircumcirclePrecisely(points);

                // Add tri as the first triad, and the three points to the convex hull
                triads.Add(tri);
                hull.Add(new HullVertex(points, tri.a));
                hull.Add(new HullVertex(points, tri.b));
                hull.Add(new HullVertex(points, tri.c));

                // Sort the remainder according to their distance from its centroid
                // Re-measure the points' distances from the centre of the circumcircle
                Vertex centre = new Vertex(circumCentreX, circumCentreY);
                for (int k = 3; k < nump; k++)
                    distance2ToCentre[k] = points[sortedIndices[k]].distance2To(centre);

                // Sort the _other_ points in order of distance to circumcentre
                Array.Sort(distance2ToCentre, sortedIndices, 3, nump - 3);

                // Add new points into hull (removing obscured ones from the chain)
                // and creating triangles....
                int numt = 0;
                for (int k = 3; k < nump; k++)
                {
                    int pointsIndex = sortedIndices[k];
                    HullVertex ptx = new HullVertex(points, pointsIndex);

                    double dx = ptx.x - hull[0].x, dy = ptx.y - hull[0].y;  // outwards pointing from hull[0] to pt.

                    int numh = hull.Count, numh_old = numh;
                    List<int> pidx = new List<int>(), tridx = new List<int>();
                    int hidx;  // new hull point location within hull.....

                    if (hull.EdgeVisibleFrom(0, dx, dy))
                    {
                        // starting with a visible hull facet !!!
                        int e2 = numh;
                        hidx = 0;

                        // check to see if segment numh is also visible
                        if (hull.EdgeVisibleFrom(numh - 1, dx, dy))
                        {
                            // visible.
                            pidx.Add(hull[numh - 1].pointsIndex);
                            tridx.Add(hull[numh - 1].triadIndex);

                            for (int h = 0; h < numh - 1; h++)
                            {
                                // if segment h is visible delete h
                                pidx.Add(hull[h].pointsIndex);
                                tridx.Add(hull[h].triadIndex);
                                if (hull.EdgeVisibleFrom(h, ptx))
                                {
                                    hull.RemoveAt(h);
                                    h--;
                                    numh--;
                                }
                                else
                                {
                                    // quit on invisibility
                                    hull.Insert(0, ptx);
                                    numh++;
                                    break;
                                }
                            }
                            // look backwards through the hull structure
                            for (int h = numh - 2; h > 0; h--)
                            {
                                // if segment h is visible delete h + 1
                                if (hull.EdgeVisibleFrom(h, ptx))
                                {
                                    pidx.Insert(0, hull[h].pointsIndex);
                                    tridx.Insert(0, hull[h].triadIndex);
                                    hull.RemoveAt(h + 1);  // erase end of chain
                                }
                                else
                                    break; // quit on invisibility
                            }
                        }
                        else
                        {
                            hidx = 1;  // keep pt hull[0]
                            tridx.Add(hull[0].triadIndex);
                            pidx.Add(hull[0].pointsIndex);

                            for (int h = 1; h < numh; h++)
                            {
                                // if segment h is visible delete h  
                                pidx.Add(hull[h].pointsIndex);
                                tridx.Add(hull[h].triadIndex);
                                if (hull.EdgeVisibleFrom(h, ptx))
                                {                     // visible
                                    hull.RemoveAt(h);
                                    h--;
                                    numh--;
                                }
                                else
                                {
                                    // quit on invisibility
                                    hull.Insert(h, ptx);
                                    break;
                                }
                            }
                        }
                    }
                    else
                    {
                        int e1 = -1, e2 = numh;
                        for (int h = 1; h < numh; h++)
                        {
                            if (hull.EdgeVisibleFrom(h, ptx))
                            {
                                if (e1 < 0)
                                    e1 = h;  // first visible
                            }
                            else
                            {
                                if (e1 > 0)
                                {
                                    // first invisible segment.
                                    e2 = h;
                                    break;
                                }
                            }
                        }

                        // triangle pidx starts at e1 and ends at e2 (inclusive).	
                        if (e2 < numh)
                        {
                            for (int e = e1; e <= e2; e++)
                            {
                                pidx.Add(hull[e].pointsIndex);
                                tridx.Add(hull[e].triadIndex);
                            }
                        }
                        else
                        {
                            for (int e = e1; e < e2; e++)
                            {
                                pidx.Add(hull[e].pointsIndex);
                                tridx.Add(hull[e].triadIndex);   // there are only n-1 triangles from n hull pts.
                            }
                            pidx.Add(hull[0].pointsIndex);
                        }

                        // erase elements e1+1 : e2-1 inclusive.
                        if (e1 < e2 - 1)
                            hull.RemoveRange(e1 + 1, e2 - e1 - 1);

                        // insert ptx at location e1+1.
                        hull.Insert(e1 + 1, ptx);
                        hidx = e1 + 1;
                    }

                    // If we're only computing the hull, we're done with this point
                    if (hullOnly)
                        continue;

                    int a = pointsIndex, T0;

                    int npx = pidx.Count - 1;
                    numt = triads.Count;
                    T0 = numt;

                    for (int p = 0; p < npx; p++)
                    {
                        Triad trx = new Triad(a, pidx[p], pidx[p + 1]);
                        trx.FindCircumcirclePrecisely(points);

                        trx.bc = tridx[p];
                        if (p > 0)
                            trx.ab = numt - 1;
                        trx.ac = numt + 1;

                        // index back into the triads.
                        Triad txx = triads[tridx[p]];
                        if ((trx.b == txx.a && trx.c == txx.b) | (trx.b == txx.b && trx.c == txx.a))
                            txx.ab = numt;
                        else if ((trx.b == txx.a && trx.c == txx.c) | (trx.b == txx.c && trx.c == txx.a))
                            txx.ac = numt;
                        else if ((trx.b == txx.b && trx.c == txx.c) | (trx.b == txx.c && trx.c == txx.b))
                            txx.bc = numt;

                        triads.Add(trx);
                        numt++;
                    }
                    // Last edge is on the outside
                    triads[numt - 1].ac = -1;

                    hull[hidx].triadIndex = numt - 1;
                    if (hidx > 0)
                        hull[hidx - 1].triadIndex = T0;
                    else
                    {
                        numh = hull.Count;
                        hull[numh - 1].triadIndex = T0;
                    }
                }
            }

            /// <summary>
            /// Return the convex hull of the supplied points,
            /// Don't check for duplicate points
            /// </summary>
            /// <param name="points">List of 2D vertices</param>
            /// <returns></returns>
            public List<Vertex> ConvexHull(List<Vertex> points)
            {
                return ConvexHull(points, false);
            }

            /// <summary>
            /// Return the convex hull of the supplied points,
            /// Optionally check for duplicate points
            /// </summary>
            /// <param name="points">List of 2D vertices</param>
            /// <param name="rejectDuplicatePoints">Whether to omit duplicated points</param>
            /// <returns></returns>
            public List<Vertex> ConvexHull(List<Vertex> points, bool rejectDuplicatePoints)
            {
                Hull hull = new Hull();
                List<Triad> triads = new List<Triad>();

                Analyse(points, hull, triads, rejectDuplicatePoints, true);

                List<Vertex> hullVertices = new List<Vertex>();
                foreach (HullVertex hv in hull)
                    hullVertices.Add(new Vertex(hv.x, hv.y));

                return hullVertices;
            }

            /// <summary>
            /// Return the Delaunay triangulation of the supplied points
            /// Don't check for duplicate points
            /// </summary>
            /// <param name="points">List of 2D vertices</param>
            /// <returns>Triads specifying the triangulation</returns>
            public List<Triad> Triangulation(List<Vertex> points)
            {
                return Triangulation(points, false);
            }

            /// <summary>
            /// Return the Delaunay triangulation of the supplied points
            /// Optionally check for duplicate points
            /// </summary>
            /// <param name="points">List of 2D vertices</param>
            /// <param name="rejectDuplicatePoints">Whether to omit duplicated points</param>
            /// <returns></returns>
            public List<Triad> Triangulation(List<Vertex> points, bool rejectDuplicatePoints)
            {
                List<Triad> triads = new List<Triad>();
                Hull hull = new Hull();

                Analyse(points, hull, triads, rejectDuplicatePoints, false);

                // Now, need to flip any pairs of adjacent triangles not satisfying
                // the Delaunay criterion
                int numt = triads.Count;
                bool[] idsA = new bool[numt];
                bool[] idsB = new bool[numt];

                // We maintain a "list" of the triangles we've flipped in order to propogate any
                // consequent changes
                // When the number of changes is large, this is best maintained as a vector of bools
                // When the number becomes small, it's best maintained as a set
                // We switch between these regimes as the number flipped decreases
                //
                // the iteration cycle limit is included to prevent degenerate cases 'oscillating'
                // and the algorithm failing to stop.
                int flipped = FlipTriangles(triads, idsA);

                int iterations = 1;
                while (flipped > (int)(fraction * (double)numt) && iterations < 1000)
                {
                    if ((iterations & 1) == 1)
                        flipped = FlipTriangles(triads, idsA, idsB);
                    else
                        flipped = FlipTriangles(triads, idsB, idsA);

                    iterations++;
                }

                Set<int> idSetA = new Set<int>(), idSetB = new Set<int>();
                flipped = FlipTriangles(triads,
                    ((iterations & 1) == 1) ? idsA : idsB, idSetA);

                iterations = 1;
                while (flipped > 0 && iterations < 2000)
                {
                    if ((iterations & 1) == 1)
                        flipped = FlipTriangles(triads, idSetA, idSetB);
                    else
                        flipped = FlipTriangles(triads, idSetB, idSetA);

                    iterations++;
                }

                return triads;
            }

            public double fraction = 0.3f;

            /// <summary>
            /// Test the triad against its 3 neighbours and flip it with any neighbour whose opposite point
            /// is inside the circumcircle of the triad
            /// </summary>
            /// <param name="triads">The triads</param>
            /// <param name="triadIndexToTest">The index of the triad to test</param>
            /// <param name="triadIndexFlipped">Index of adjacent triangle it was flipped with (if any)</param>
            /// <returns>true iff the triad was flipped with any of its neighbours</returns>
            bool FlipTriangle(List<Triad> triads, int triadIndexToTest, out int triadIndexFlipped)
            {
                int oppositeVertex = 0, edge1, edge2, edge3 = 0, edge4 = 0;
                triadIndexFlipped = 0;

                Triad tri = triads[triadIndexToTest];
                // test all 3 neighbours of tri 

                if (tri.bc >= 0)
                {
                    triadIndexFlipped = tri.bc;
                    Triad t2 = triads[triadIndexFlipped];
                    // find relative orientation (shared limb).
                    t2.FindAdjacency(tri.b, triadIndexToTest, out oppositeVertex, out edge3, out edge4);
                    if (tri.InsideCircumcircle(points[oppositeVertex]))
                    {  // not valid in the Delaunay sense.
                        edge1 = tri.ab;
                        edge2 = tri.ac;
                        if (edge1 != edge3 && edge2 != edge4)
                        {
                            int tria = tri.a, trib = tri.b, tric = tri.c;
                            tri.Initialize(tria, trib, oppositeVertex, edge1, edge3, triadIndexFlipped, points);
                            t2.Initialize(tria, tric, oppositeVertex, edge2, edge4, triadIndexToTest, points);

                            // change knock on triangle labels.
                            if (edge3 >= 0)
                                triads[edge3].ChangeAdjacentIndex(triadIndexFlipped, triadIndexToTest);
                            if (edge2 >= 0)
                                triads[edge2].ChangeAdjacentIndex(triadIndexToTest, triadIndexFlipped);
                            return true;
                        }
                    }
                }


                if (tri.ab >= 0)
                {
                    triadIndexFlipped = tri.ab;
                    Triad t2 = triads[triadIndexFlipped];
                    // find relative orientation (shared limb).
                    t2.FindAdjacency(tri.a, triadIndexToTest, out oppositeVertex, out edge3, out edge4);
                    if (tri.InsideCircumcircle(points[oppositeVertex]))
                    {  // not valid in the Delaunay sense.
                        edge1 = tri.ac;
                        edge2 = tri.bc;
                        if (edge1 != edge3 && edge2 != edge4)
                        {
                            int tria = tri.a, trib = tri.b, tric = tri.c;
                            tri.Initialize(tric, tria, oppositeVertex, edge1, edge3, triadIndexFlipped, points);
                            t2.Initialize(tric, trib, oppositeVertex, edge2, edge4, triadIndexToTest, points);

                            // change knock on triangle labels.
                            if (edge3 >= 0)
                                triads[edge3].ChangeAdjacentIndex(triadIndexFlipped, triadIndexToTest);
                            if (edge2 >= 0)
                                triads[edge2].ChangeAdjacentIndex(triadIndexToTest, triadIndexFlipped);
                            return true;
                        }
                    }
                }

                if (tri.ac >= 0)
                {
                    triadIndexFlipped = tri.ac;
                    Triad t2 = triads[triadIndexFlipped];
                    // find relative orientation (shared limb).
                    t2.FindAdjacency(tri.a, triadIndexToTest, out oppositeVertex, out edge3, out edge4);
                    if (tri.InsideCircumcircle(points[oppositeVertex]))
                    {  // not valid in the Delaunay sense.
                        edge1 = tri.ab;   // .ac shared limb
                        edge2 = tri.bc;
                        if (edge1 != edge3 && edge2 != edge4)
                        {
                            int tria = tri.a, trib = tri.b, tric = tri.c;
                            tri.Initialize(trib, tria, oppositeVertex, edge1, edge3, triadIndexFlipped, points);
                            t2.Initialize(trib, tric, oppositeVertex, edge2, edge4, triadIndexToTest, points);

                            // change knock on triangle labels.
                            if (edge3 >= 0)
                                triads[edge3].ChangeAdjacentIndex(triadIndexFlipped, triadIndexToTest);
                            if (edge2 >= 0)
                                triads[edge2].ChangeAdjacentIndex(triadIndexToTest, triadIndexFlipped);
                            return true;
                        }
                    }
                }

                return false;
            }

            /// <summary>
            /// Flip triangles that do not satisfy the Delaunay condition
            /// </summary>
            private int FlipTriangles(List<Triad> triads, bool[] idsFlipped)
            {
                int numt = (int)triads.Count;
                Array.Clear(idsFlipped, 0, numt);

                int flipped = 0;
                for (int t = 0; t < numt; t++)
                {
                    int t2;
                    if (FlipTriangle(triads, t, out t2))
                    {
                        flipped += 2;
                        idsFlipped[t] = true;
                        idsFlipped[t2] = true;

                    }
                }

                return flipped;
            }

            private int FlipTriangles(List<Triad> triads, bool[] idsToTest, bool[] idsFlipped)
            {
                int numt = (int)triads.Count;
                Array.Clear(idsFlipped, 0, numt);

                int flipped = 0;
                for (int t = 0; t < numt; t++)
                {
                    if (idsToTest[t])
                    {
                        int t2;
                        if (FlipTriangle(triads, t, out t2))
                        {
                            flipped += 2;
                            idsFlipped[t] = true;
                            idsFlipped[t2] = true;
                        }
                    }
                }

                return flipped;
            }

            private int FlipTriangles(List<Triad> triads, bool[] idsToTest, Set<int> idsFlipped)
            {
                int numt = (int)triads.Count;
                idsFlipped.Clear();

                int flipped = 0;
                for (int t = 0; t < numt; t++)
                {
                    if (idsToTest[t])
                    {
                        int t2;
                        if (FlipTriangle(triads, t, out t2))
                        {
                            flipped += 2;
                            idsFlipped.Add(t);
                            idsFlipped.Add(t2);
                        }
                    }
                }

                return flipped;
            }

            private int FlipTriangles(List<Triad> triads, Set<int> idsToTest, Set<int> idsFlipped)
            {
                int flipped = 0;
                idsFlipped.Clear();

                foreach (int t in idsToTest)
                {
                    int t2;
                    if (FlipTriangle(triads, t, out t2))
                    {
                        flipped += 2;
                        idsFlipped.Add(t);
                        idsFlipped.Add(t2);
                    }
                }

                return flipped;
            }

            #region Debug verification routines: verify that triad adjacency and indeces are set correctly
#if DEBUG
            private void VerifyHullContains(Hull hull, int idA, int idB)
            {
                if (
                    ((hull[0].pointsIndex == idA) && (hull[hull.Count - 1].pointsIndex == idB)) ||
                    ((hull[0].pointsIndex == idB) && (hull[hull.Count - 1].pointsIndex == idA)))
                    return;

                for (int h = 0; h < hull.Count - 1; h++)
                {
                    if (hull[h].pointsIndex == idA)
                    {
                        Debug.Assert(hull[h + 1].pointsIndex == idB);
                        return;
                    }
                    else if (hull[h].pointsIndex == idB)
                    {
                        Debug.Assert(hull[h + 1].pointsIndex == idA);
                        return;
                    }
                }

            }

            private void VerifyTriadContains(Triad tri, int nbourTriad, int idA, int idB)
            {
                if (tri.ab == nbourTriad)
                {
                    Debug.Assert(
                        ((tri.a == idA) && (tri.b == idB)) ||
                        ((tri.b == idA) && (tri.a == idB)));
                }
                else if (tri.ac == nbourTriad)
                {
                    Debug.Assert(
                        ((tri.a == idA) && (tri.c == idB)) ||
                        ((tri.c == idA) && (tri.a == idB)));
                }
                else if (tri.bc == nbourTriad)
                {
                    Debug.Assert(
                        ((tri.c == idA) && (tri.b == idB)) ||
                        ((tri.b == idA) && (tri.c == idB)));
                }
                else
                    Debug.Assert(false);
            }

            private void VerifyTriads(List<Triad> triads, Hull hull)
            {
                for (int t = 0; t < triads.Count; t++)
                {
                    if (t == 17840)
                        t = t + 0;

                    Triad tri = triads[t];
                    if (tri.ac == -1)
                        VerifyHullContains(hull, tri.a, tri.c);
                    else
                        VerifyTriadContains(triads[tri.ac], t, tri.a, tri.c);

                    if (tri.ab == -1)
                        VerifyHullContains(hull, tri.a, tri.b);
                    else
                        VerifyTriadContains(triads[tri.ab], t, tri.a, tri.b);

                    if (tri.bc == -1)
                        VerifyHullContains(hull, tri.b, tri.c);
                    else
                        VerifyTriadContains(triads[tri.bc], t, tri.b, tri.c);

                }
            }

            private void WriteTriangles(List<Triad> triangles, string name)
            {
                using (System.IO.StreamWriter writer = new System.IO.StreamWriter(name + ".dtt"))
                {
                    writer.WriteLine(triangles.Count.ToString());
                    for (int i = 0; i < triangles.Count; i++)
                    {
                        Triad t = triangles[i];
                        writer.WriteLine(string.Format("{0}: {1} {2} {3} - {4} {5} {6}",
                            i + 1,
                            t.a, t.b, t.c,
                            t.ab + 1, t.bc + 1, t.ac + 1));
                    }
                }
            }

#endif

            #endregion
        }

        class Triad
        {
            public int a, b, c;     // номера точек массива (номера вершин)
            public int ab, bc, ac;  // adjacent edges index to neighbouring triangle.

            // Position and radius squared of circumcircle
            public double circumcircleR2, circumcircleX, circumcircleY;

            public Triad(int x, int y, int z)
            {
                a = x; b = y; c = z; ab = -1; bc = -1; ac = -1;
                circumcircleR2 = -1; //x = 0; y = 0;
            }

            public void Initialize(int a, int b, int c, int ab, int bc, int ac, List<Vertex> points)
            {
                this.a = a;
                this.b = b;
                this.c = c;
                this.ab = ab;
                this.bc = bc;
                this.ac = ac;

                FindCircumcirclePrecisely(points);
            }

            /// <summary>
            /// If current orientation is not clockwise, swap b<->c
            /// </summary>
            public void MakeClockwise(List<Vertex> points)
            {
                double centroidX = (points[a].x + points[b].x + points[c].x) / 3.0f;
                double centroidY = (points[a].y + points[b].y + points[c].y) / 3.0f;

                double dr0 = points[a].x - centroidX, dc0 = points[a].y - centroidY;
                double dx01 = points[b].x - points[a].x, dy01 = points[b].y - points[a].y;

                double df = -dx01 * dc0 + dy01 * dr0;
                if (df > 0)
                {
                    // Need to swap vertices b<->c and edges ab<->bc
                    int t = b;
                    b = c;
                    c = t;

                    t = ab;
                    ab = ac;
                    ac = t;
                }
            }

            /// <summary>
            /// Find location and radius ^2 of the circumcircle (through all 3 points)
            /// This is the most critical routine in the entire set of code.  It must
            /// be numerically stable when the points are nearly collinear.
            /// </summary>
            public bool FindCircumcirclePrecisely(List<Vertex> points)
            {
                // Use coordinates relative to point `a' of the triangle
                Vertex pa = points[a], pb = points[b], pc = points[c];

                double xba = pb.x - pa.x;
                double yba = pb.y - pa.y;
                double xca = pc.x - pa.x;
                double yca = pc.y - pa.y;

                // Squares of lengths of the edges incident to `a'
                double balength = xba * xba + yba * yba;
                double calength = xca * xca + yca * yca;

                // Calculate the denominator of the formulae. 
                double D = xba * yca - yba * xca;
                if (D == 0)
                {
                    circumcircleX = 0;
                    circumcircleY = 0;
                    circumcircleR2 = -1;
                    return false;
                }

                double denominator = 0.5 / D;

                // Calculate offset (from pa) of circumcenter
                double xC = (yca * balength - yba * calength) * denominator;
                double yC = (xba * calength - xca * balength) * denominator;

                double radius2 = xC * xC + yC * yC;
                if ((radius2 > 1e10 * balength || radius2 > 1e10 * calength))
                {
                    circumcircleX = 0;
                    circumcircleY = 0;
                    circumcircleR2 = -1;
                    return false;
                }

                circumcircleR2 = (double)radius2;
                circumcircleX = (double)(pa.x + xC);
                circumcircleY = (double)(pa.y + yC);

                return true;
            }

            /// <summary>
            /// Return true iff Vertex p is inside the circumcircle of this triangle
            /// </summary>
            public bool InsideCircumcircle(Vertex p)
            {
                double dx = circumcircleX - p.x;
                double dy = circumcircleY - p.y;
                double r2 = dx * dx + dy * dy;
                return r2 < circumcircleR2;
            }

            /// <summary>
            /// Change any adjacent triangle index that matches fromIndex, to toIndex
            /// </summary>
            public void ChangeAdjacentIndex(int fromIndex, int toIndex)
            {
                if (ab == fromIndex)
                    ab = toIndex;
                else if (bc == fromIndex)
                    bc = toIndex;
                else if (ac == fromIndex)
                    ac = toIndex;
                else
                    Debug.Assert(false);
            }

            /// <summary>
            /// Determine which edge matches the triangleIndex, then which vertex the vertexIndex
            /// Set the indices of the opposite vertex, left and right edges accordingly
            /// </summary>
            public void FindAdjacency(int vertexIndex, int triangleIndex, out int indexOpposite, out int indexLeft, out int indexRight)
            {
                if (ab == triangleIndex)
                {
                    indexOpposite = c;

                    if (vertexIndex == a)
                    {
                        indexLeft = ac;
                        indexRight = bc;
                    }
                    else
                    {
                        indexLeft = bc;
                        indexRight = ac;
                    }
                }
                else if (ac == triangleIndex)
                {
                    indexOpposite = b;

                    if (vertexIndex == a)
                    {
                        indexLeft = ab;
                        indexRight = bc;
                    }
                    else
                    {
                        indexLeft = bc;
                        indexRight = ab;
                    }
                }
                else if (bc == triangleIndex)
                {
                    indexOpposite = a;

                    if (vertexIndex == b)
                    {
                        indexLeft = ab;
                        indexRight = ac;
                    }
                    else
                    {
                        indexLeft = ac;
                        indexRight = ab;
                    }
                }
                else
                {
                    Debug.Assert(false);
                    indexOpposite = indexLeft = indexRight = 0;
                }
            }

            public override string ToString()
            {
                return string.Format("Triad vertices {0} {1} {2} ; edges {3} {4} {5}", a, b, c, ab, ac, bc);
            }
        }

        class Set<T> : IEnumerable<T>
        {
            SortedList<T, int> list;

            public Set()
            {
                list = new SortedList<T, int>();
            }

            public void Add(T k)
            {
                if (!list.ContainsKey(k))
                    list.Add(k, 0);
            }

            public int Count
            {
                get { return list.Count; }
            }

            public void DeepCopy(Set<T> other)
            {
                list.Clear();
                foreach (T k in other.list.Keys)
                    Add(k);
            }

            IEnumerator<T> IEnumerable<T>.GetEnumerator()
            {
                return list.Keys.GetEnumerator();
            }

            public void Clear()
            {
                list.Clear();
            }


            #region IEnumerable Members

            System.Collections.IEnumerator System.Collections.IEnumerable.GetEnumerator()
            {
                throw new NotImplementedException();
            }

            #endregion
        }

        /// <summary>
        /// Vertices belonging to the convex hull need to maintain a point and triad index
        /// </summary>
        internal class HullVertex : Vertex
        {
            public int pointsIndex;
            public int triadIndex;

            public HullVertex(List<Vertex> points, int pointIndex)
            {
                x = points[pointIndex].x;
                y = points[pointIndex].y;
                pointsIndex = pointIndex;
                triadIndex = 0;
            }
        }

        /// <summary>
        /// Hull represents a list of vertices in the convex hull, and keeps track of
        /// their indices (into the associated points list) and triads
        /// </summary>
        class Hull : List<HullVertex>
        {
            private int NextIndex(int index)
            {
                if (index == Count - 1)
                    return 0;
                else
                    return index + 1;
            }

            /// <summary>
            /// Return vector from the hull point at index to next point
            /// </summary>
            public void VectorToNext(int index, out double dx, out double dy)
            {
                Vertex et = this[index], en = this[NextIndex(index)];

                dx = en.x - et.x;
                dy = en.y - et.y;
            }

            /// <summary>
            /// Return whether the hull vertex at index is visible from the supplied coordinates
            /// </summary>
            public bool EdgeVisibleFrom(int index, double dx, double dy)
            {
                double idx, idy;
                VectorToNext(index, out idx, out idy);

                double crossProduct = -dy * idx + dx * idy;
                return crossProduct < 0;
            }

            /// <summary>
            /// Return whether the hull vertex at index is visible from the point
            /// </summary>
            public bool EdgeVisibleFrom(int index, Vertex point)
            {
                double idx, idy;
                VectorToNext(index, out idx, out idy);

                double dx = point.x - this[index].x;
                double dy = point.y - this[index].y;

                double crossProduct = -dy * idx + dx * idy;
                return crossProduct < 0;
            }
        }
    }

    namespace GeometryMath
    {
        /// <summary>
        /// Точка.
        /// </summary>
        internal class Point
        {
            public double x;
            public double y;
            public int i; // дополнительный параметр - индекс точки
            // Создание точки (_x, _y)
            public Point(double _x, double _y)
            {
                x = _x;
                y = _y;
            }
            // Создание точки (0, 0)
            public Point()
            {
                x = 0;
                y = 0;
            }

            public static bool operator ==(Point a, Point b)
            {
                if (a.x == b.x && a.y == b.y)
                    return true;
                return false;
            }

            public static bool operator !=(Point a, Point b)
            {
                if (a.x != b.x || a.y != b.y)
                    return true;
                return false;
            }
        }

        /// <summary>
        /// Прямая.
        /// </summary>
        internal class Line
        {
            public double a;
            public double b;
            public double c;
            // Создание прямой ax + by + c = 0
            public Line(double _a = 0, double _b = 0, double _c = 0)
            {
                a = _a;
                b = _b;
                c = _c;
            }
        }

        /// <summary>
        /// Окружность.
        /// </summary>
        internal class Circle
        {
            public Point c;     // центр
            public double r;    // радиус
            // дополнительные данные:
            public double alpha;
            // Создание окружности с центром (x, y) и радиусом _r
            public Circle (double x, double y, double _r)
            {
               c = new Point(x, y);
               r = _r;
            }
            // Создание окружности с центром p и радиусом _r
            public Circle (Point p, double _r)
            {
                c = new Point(p.x, p.y);
                r = _r;
            }
            public Circle ()
            {
                c = new Point(0, 0);
                r = 0;
            }
        }

        /// <summary>
        /// Отрезок.
        /// </summary>
        internal class Segment
        {
            public Point a;
            public Point b;

            public Segment()
            {
                a = new Point(0, 0);
                b = new Point(0, 0);
            }
            public Segment(Point _a, Point _b)
            {
                a = new Point(_a.x, _a.y);
                b = new Point(_b.x, _b.y);
            }
        }


        internal class Triangle
        {
            public Point a;
            public Point b;
            public Point c;

            public Triangle()
            {
                a = new Point(0, 0);
                b = new Point(0, 0);
                c = new Point(0, 0);
            }
            public Triangle(Point _a, Point _b, Point _c)
            {
                a = new Point(_a.x, _a.y);
                b = new Point(_b.x, _b.y);
                c = new Point(_c.x, _c.y);
            }
        }

        /// <summary>
        /// Функции работы с точками.
        /// </summary>
        internal static class funcPoint
        {
            /// <summary>
            /// Расстояние между двумя точками.
            /// </summary>
            internal static double dist(Point p1, Point p2)
            {
                return Math.Sqrt((p1.x - p2.x) * (p1.x - p2.x) + (p1.y - p2.y) * (p1.y - p2.y));
            }

            /// <summary>
            /// Принадлежноть точки прямоугольной области
            /// </summary>
            internal static bool point_in_box(Point t, Point p1, Point p2)
            {
                const double eps = 1e-8;
                return (Math.Abs(t.x - Math.Min(p1.x, p2.x)) <= eps || Math.Min(p1.x, p2.x) <= t.x) &&
                        (Math.Abs(Math.Max(p1.x, p2.x) - t.x) <= eps || Math.Max(p1.x, p2.x) >= t.x) &&
                        (Math.Abs(t.y - Math.Min(p1.y, p2.y)) <= eps || Math.Min(p1.y, p2.y) <= t.y) &&
                        (Math.Abs(Math.Max(p1.y, p2.y) - t.y) <= eps || Math.Max(p1.y, p2.y) >= t.y);
            }

            /// <summary>
            /// Деление отрезка в заданном соотношении.
            /// </summary>
            /// <returns>Точка разрыва.</returns>
            internal static Point part_segment(Point p1, Point p2, double m, double n)
            {
                Point t = new Point();
                t.x = (p1.x * n + p2.x * m) / (m + n);
                t.y = (p1.y * n + p2.y * m) / (m + n);
                return t;
            }
        }

        /// <summary>
        /// Функции работы с прямыми.
        /// </summary>
        internal static class funcLine
        {
            /// <summary>
            /// Уравнение прямой, проходящей через 2 точки.
            /// </summary>
            internal static Line toline(Point p1, Point p2)
            {
                double a = p2.y - p1.y;
                double b = p1.x - p2.x;

                return new Line(a, b, -a * p1.x - b * p1.y);
            }

            /// <summary>
            /// Ориентация точки относительно прямой.
            /// </summary>
            internal static int point_in_line(Line l, Point p)
            {
                const double eps = 1e-8;
                double s = l.a * p.x + l.b * p.y + l.c;
                return s < -eps ? -1 : s > eps ? 1 : 0;
            }

            /// <summary>
            /// Ориентация точки относительно прямой. (Точки). -1: слева, 0: на прямой, 1: справа.
            /// </summary>
            internal static int point_in_line_ADVANCED(Point line_point1, Point line_point2, Point testPoint)
            {
                double tmp = (line_point2.x - line_point1.x) * (testPoint.y - line_point1.y) -
                          (line_point2.y - line_point1.y) * (testPoint.x - line_point1.x);

                if (tmp > 0)
                    return 1;
                else
                    if (tmp < 0)
                        return -1;
                    else
                        return 0;
            }

            /// <summary>
            /// Параллельность прямых.
            /// </summary>
            internal static bool is_parallel_line(Line l1, Line l2)
            {
                const double eps = 1e-8;
                return Math.Abs(l1.a * l2.b - l2.a * l1.b) <= eps;
            }

            /// <summary>
            /// Совпадение прямых.
            /// </summary>
            internal static bool is_equal_line(Line l1, Line l2)
            {
                const double eps = 1e-8;
                return Math.Abs(l1.a * l2.b - l2.a * l1.b) <= eps &&
                        Math.Abs(l1.a * l2.c - l2.a * l1.c) <= eps &&
                        Math.Abs(l1.b * l2.c - l2.b * l1.c) <= eps;
            }

            /// <summary>
            /// Пересекаются ли прямые?
            /// </summary>
            internal static bool is_cross_line(Point p1, Point p2, Point p3, Point p4)
            {
                Line l1 = funcLine.toline(p1, p2);
                Line l2 = funcLine.toline(p3, p4);

                if (is_equal_line(l1, l2))
                    return false;
                if (is_parallel_line(l1, l2))
                    return false;
                return true;
            }

            /// <summary>
            /// Точка пересечения прямых.
            /// </summary>
            /// <returns>
            /// Если прямые параллельные, функция возвращает 0.
            /// Если прямые совпадают, функция возвращает 2.
            /// Иначе, прямые пересекаются и функция возвращет 1.
            /// </returns>
            internal static int cross_line(Line l1, Line l2, ref Point p)
            {
                if (is_equal_line (l1, l2)) 
                    return 2;
                if (is_parallel_line (l1, l2)) 
                    return 0;

                p.x = (l2.b * l1.c - l1.b * l2.c) / (l2.a * l1.b - l1.a * l2.b);
                p.y = (l1.b != 0 ? (- l1.c - l1.a * p.x) / l1.b : (- l2.c - l2.a * p.x) / l2.b);

                return 1;
            }

            /// <summary>
            /// Перпендикуляр к прямой.
            /// </summary>
            internal static Line perp_line(Line l, Point p)
            {
                return new Line(l.b, -l.a, -l.b * p.x + l.a * p.y);
            }

            /// <summary>
            /// Pасстояние от точки до прямой.
            /// </summary>
            internal static double dist_point_to_line(Point p, Line l)
            {
                return Math.Abs((l.a * p.x + l.b * p.y + l.c) / Math.Sqrt(l.a * l.a + l.b * l.b));
            }

            /// <summary>
            /// Pасстояние между параллельными прямыми
            /// </summary>
            internal static double dist_between_line(Line l1, Line l2)
            {
                return Math.Abs(l1.c - l2.c) / Math.Sqrt(l1.a * l1.a + l2.b * l2.b);
            }
        }

        /// <summary>
        /// Функции работы с отрезками.
        /// </summary>
        internal static class funcSegment
        {
            /// <summary>
            /// Идентичны ли отрезки?
            /// </summary>
            internal static bool is_identical(Segment s1, Segment s2)
            {
                if ((s1.a == s2.a && s1.b == s2.b) || (s1.a == s2.b && s1.b == s2.a))
                    return true;
                return false;
            }

            /// <summary>
            /// Принадлежность точки отрезку.
            /// </summary>
            internal static bool point_in_segment(Point t, Point p1, Point p2)
            {
                const double eps = 1e-8;
                double a = p2.y - p1.y;
                double b = p1.x - p2.x;
                double c = -a * p1.x - b * p1.y;
                if (Math.Abs(a * t.x + b * t.y + c) > eps) 
                    return false;

                return funcPoint.point_in_box(t, p1, p2);
            }

            /// <summary>
            /// Пересекаются ли отрезки?
            /// </summary>
            internal static bool is_cross_segment(Point p1, Point p2, Point p3, Point p4)
            {
                const double eps = 1e-8;
                Line l1 = funcLine.toline(p1, p2);
                Line l2 = funcLine.toline(p3, p4);
                int sign1 = funcLine.point_in_line(l1, p3) * funcLine.point_in_line(l1, p4);
                int sign2 = funcLine.point_in_line(l2, p1) * funcLine.point_in_line(l2, p2);

                if (Math.Abs(sign1) <= eps && Math.Abs(sign2) <= eps)
                    return funcPoint.point_in_box(p1, p3, p4) || funcPoint.point_in_box(p2, p3, p4) ||
                           funcPoint.point_in_box(p3, p1, p2) || funcPoint.point_in_box(p4, p1, p2);
                return sign1 <= eps && sign2 <= eps;
            }

            /// <summary>
            /// Пересекаются ли отрезки? Крайные точки выколоты.
            /// </summary>
            internal static bool is_cross_segment_ADVANCED(Point p1, Point p2, Point p3, Point p4)
            {
                if (p1 == p3 || p1 == p4 || p2 == p3 || p2 == p4)
                    return false;

                const double eps = 1e-8;
                Line l1 = funcLine.toline(p1, p2);
                Line l2 = funcLine.toline(p3, p4);
                int sign1 = funcLine.point_in_line(l1, p3) * funcLine.point_in_line(l1, p4);
                int sign2 = funcLine.point_in_line(l2, p1) * funcLine.point_in_line(l2, p2);

                if (Math.Abs(sign1) <= eps && Math.Abs(sign2) <= eps)
                    return funcPoint.point_in_box(p1, p3, p4) || funcPoint.point_in_box(p2, p3, p4) ||
                           funcPoint.point_in_box(p3, p1, p2) || funcPoint.point_in_box(p4, p1, p2);
                return sign1 <= eps && sign2 <= eps;
            }

            /// <summary>
            /// Точка пересечения отрезков.
            /// </summary>
            internal static bool cross_segment(Point p1, Point p2, Point p3, Point p4, ref Point t)
            {
                // Строим прямые проходящие через эти отрезки и пересекаем их
                Line l1 = funcLine.toline(p1, p2);
                Line l2 = funcLine.toline(p3, p4);

                int flag = funcLine.cross_line(l1, l2, ref t);
                if (flag == 0) return false;

                // Если прямые совпадают, проверяем каждый конец отрезка на принадлежность другому отрезку
                if (flag == 2)
                {
                    if (funcPoint.point_in_box(p1, p3, p4)) { t = p1; return true; }
                    if (funcPoint.point_in_box(p2, p3, p4)) { t = p2; return true; }
                    if (funcPoint.point_in_box(p3, p1, p2)) { t = p3; return true; }
                    if (funcPoint.point_in_box(p4, p1, p2)) { t = p4; return true; }
                    return false;
                }
                // Если прямые пересекаются, проверяем принадлежит ли точка пересечения обоим отрезкам
                return funcPoint.point_in_box(t, p1, p2) && funcPoint.point_in_box(t, p3, p4);
            }

            /// <summary>
            /// Пересекает ли отрезок прямую?
            /// </summary>
            /// <param name="p1">Начальная точка отрезка.</param>
            /// <param name="p2">Конечная точка отрезка.</param>
            /// <param name="p3">Начальная точка прямой.</param>
            /// <param name="p4">Конечная точка прямой</param>
            internal static bool is_cross_segment_line_MODIFY(Point p1, Point p2, Point p3, Point p4)
            {
                Line t = funcLine.toline(p1, p2);
                Line l = funcLine.toline(p3, p4);

                // Содержит точку пересечения, если такая существует.
                Point temp = new Point();

                int flag = funcLine.cross_line(l, t, ref temp);
                if (flag == 0) return false;
                if (flag == 2) return false;

                if (funcPoint.point_in_box(temp, p1, p2)) return true;
                return false;
            }
        }

        /// <summary>
        /// Функции работы с треугольниками.
        /// </summary>
        internal static class funcTriangle
        {
            /// <summary>
            /// Высота угла.
            /// </summary>
            internal static Line height_line(Point a, Point b, Point c)
            {
                return funcLine.perp_line(funcLine.toline(a, c), b);
            }

            /// <summary>
            /// Медиана угла.
            /// </summary>
            internal static Line median_line(Point a, Point b, Point c)
            {
                return funcLine.toline(b, new Point((a.x + c.x) / 2, (a.y + c.y) / 2));
            }

            /// <summary>
            /// Биссектриса угла.
            /// </summary>
            internal static Line bisector_line(Point a, Point b, Point c)
            {
                double ab = funcPoint.dist(a, b);
                double bc = funcPoint.dist(b, c);
                Point tab = funcPoint.part_segment(b, a, bc, ab);
                Point tbc = funcPoint.part_segment(b, c, ab, bc);
                Point p = funcPoint.part_segment(tab, tbc, 1, 1);
                return funcLine.toline(b, p);
            }

            /// <summary>
            /// Вписанная в треугольник окружность.
            /// </summary>
            internal static Circle entered_circle(Point a, Point b, Point c)
            {
                Line ba = bisector_line(b, a, c);
                Line bb = bisector_line(a, b, c);

                Circle res = new Circle();
                int flag = funcLine.cross_line(ba, bb, ref res.c);
                res.r = funcLine.dist_point_to_line(res.c, funcLine.toline(a, b));
                return res;
            }

            /// <summary>
            /// Описанная около треугольника окружность.
            /// </summary>
            internal static Circle described_circle(Point a, Point b, Point c)
            {
                Point tab = funcPoint.part_segment(a, b, 1, 1);
                Point tbc = funcPoint.part_segment(b, c, 1, 1);
                Line pab = funcLine.perp_line(funcLine.toline(a, b), tab);
                Line pbc = funcLine.perp_line(funcLine.toline(b, c), tbc);

                Circle res = new Circle();
                int flag = funcLine.cross_line(pab, pbc, ref res.c);
                res.r = funcPoint.dist(res.c, a);
                return res;
            }

            /// <summary>
            /// Принадлежит ли точка треугольнику?
            /// </summary>
            internal static bool point_in_triangle(Point a, Point b, Point c, Point p)
            {
                if (funcSegment.is_cross_segment_line_MODIFY(a, p, b, c) == true) return false;
                if (funcSegment.is_cross_segment_line_MODIFY(b, p, a, c) == true) return false;
                if (funcSegment.is_cross_segment_line_MODIFY(c, p, a, b) == true) return false;
                return true;
            }

            /// <summary>
            /// Лежит ли первый треугольник внутри второго?
            /// </summary>
            internal static bool is_triangle_in_triangle(double[][] smallT, double[][] bigT)
            {
                Point smallA = new Point(smallT[0][0], smallT[0][1]);
                Point smallB = new Point(smallT[1][0], smallT[1][1]);
                Point smallC = new Point(smallT[2][0], smallT[2][1]);
                //
                Point bigA = new Point(bigT[0][0], bigT[0][1]);
                Point bigB = new Point(bigT[1][0], bigT[1][1]);
                Point bigC = new Point(bigT[2][0], bigT[2][1]);
                //
                Circle smallCircleIn = entered_circle(smallA, smallB, smallC);
                Point smallCenter = new Point(smallCircleIn.c.x, smallCircleIn.c.y);
                //
                if (point_in_triangle(bigA, bigB, bigC, smallCenter) == true)
                    return true;
                return false;
            }

            /// <summary>
            /// Лежит ли первый треугольник внутри второго (Списки)?
            /// </summary>
            internal static List<List<bool>> is_triangle_in_triangle_ADVANCED(List<double[][]> smallT, List<double[][]> bigT)
            {
                List<List<bool>> resultList = new List<List<bool>>();
                //
                for (int i = 0; i < bigT.Count; i++)
                {
                    List<bool> boolList = new List<bool>();
                    //
                    for (int j = 0; j < smallT.Count; j++)
                    {
                        if (is_triangle_in_triangle(smallT[j], bigT[i]))
                            boolList.Add(true);
                        else
                            boolList.Add(false);
                    }
                    //
                    resultList.Add(boolList);
                }
                //
                return resultList;
            }

            /// <summary>
            /// Пересекает ли отрезок ребра треугольника?
            /// </summary>
            internal static bool triangle_intersect(Triangle triangle, Point lpA, Point lpB)
            {
                if (funcSegment.is_cross_segment_ADVANCED(triangle.a, triangle.b, lpA, lpB)) return true;
                if (funcSegment.is_cross_segment_ADVANCED(triangle.b, triangle.c, lpA, lpB)) return true;
                if (funcSegment.is_cross_segment_ADVANCED(triangle.a, triangle.c, lpA, lpB)) return true;
                return false;
            }
        }
    }

    // Интерфейсы для прикладного кода.
    namespace Actions
    {
        public static class Triangulation
        {
            /// <summary>
            /// (Триангуляция) построение набора треугольников из набора вершин.
            /// </summary>
            /// <param name="triangles">
            /// List int[3 вершины][2 координаты x и y] Большие треугольники (структурные).
            /// </param>
            /// <returns>
            /// List int[3 вершины][2 координаты x и y] Маленькие треугольники (тесселяция).
            /// </returns>
            //public static void triangulate(List<double[][]> triangles)
            public static List<double[][]> triangulate(List<double[][]> triangles)
            {
                int i, j, k;
                // 1) Создать список ребер:
                i = 0;
                List<GeometryMath.Segment> structEdges = new List<GeometryMath.Segment>();
                for (i = 0; i < triangles.Count; i++)
                {
                    structEdges.Add(new GeometryMath.Segment(new GeometryMath.Point(triangles[i][0][0], triangles[i][0][1]), new GeometryMath.Point(triangles[i][1][0], triangles[i][1][1])));
                    structEdges.Add(new GeometryMath.Segment(new GeometryMath.Point(triangles[i][1][0], triangles[i][1][1]), new GeometryMath.Point(triangles[i][2][0], triangles[i][2][1])));
                    structEdges.Add(new GeometryMath.Segment(new GeometryMath.Point(triangles[i][0][0], triangles[i][0][1]), new GeometryMath.Point(triangles[i][2][0], triangles[i][2][1])));
                }

                // Округлим координаты до 2х знаков:
                i = 0; j = 0;
                for (i = 0; i < structEdges.Count; i++)
                {
                    structEdges[i].a.x = Math.Round(structEdges[i].a.x, 2);
                    structEdges[i].a.y = Math.Round(structEdges[i].a.y, 2);
                    structEdges[i].b.x = Math.Round(structEdges[i].b.x, 2);
                    structEdges[i].b.y = Math.Round(structEdges[i].b.y, 2);
                }

                // 2) Раздробить ребра в точках пересечения.
                i = -1; j = -1; // чтоб начать с 0;
                while (i < structEdges.Count - 1)
                {
                    i++;
                    while (j < structEdges.Count - 1)
                    {
                        j++;
                        if (i == j) continue;
                        // Если ребра пересекаются - создать четыре новых ребра, добавить в конец, а старые удалить.
                        if (GeometryMath.funcSegment.is_cross_segment_ADVANCED(structEdges[i].a, structEdges[i].b, structEdges[j].a, structEdges[j].b) == true)
                        {
                            // Точка пересечения
                            GeometryMath.Point p = new GeometryMath.Point();
                            GeometryMath.funcSegment.cross_segment(structEdges[i].a, structEdges[i].b, structEdges[j].a, structEdges[j].b, ref p);
                            //
                            p.x = Math.Round(p.x, 2);
                            p.y = Math.Round(p.y, 2);
                            //
                            // Создадим 4 новых ребра [ненулевые]
                            if (structEdges[i].a != p)
                            {
                                //structEdges.Add(new GeometryMath.Segment(structEdges[i].a, p));
                                structEdges.Add(
                                    new GeometryMath.Segment(
                                        new GeometryMath.Point(
                                            Math.Round(structEdges[i].a.x, 2), Math.Round(structEdges[i].a.y, 2)),
                                        new GeometryMath.Point(
                                            Math.Round(p.x, 2), Math.Round(p.y, 2))));
                            }
                            if (p != structEdges[i].b)
                            {
                                //structEdges.Add(new GeometryMath.Segment(p, structEdges[i].b));
                                structEdges.Add(
                                    new GeometryMath.Segment(
                                        new GeometryMath.Point(
                                            Math.Round(p.x, 2), Math.Round(p.y, 2)),
                                        new GeometryMath.Point(
                                            Math.Round(structEdges[i].b.x, 2), Math.Round(structEdges[i].b.y, 2))));
                            }
                            if (structEdges[j].a != p)
                            {
                                //structEdges.Add(new GeometryMath.Segment(structEdges[j].a, p));
                                structEdges.Add(
                                    new GeometryMath.Segment(
                                        new GeometryMath.Point(
                                            Math.Round(structEdges[j].a.x, 2), Math.Round(structEdges[j].a.y, 2)),
                                        new GeometryMath.Point(
                                            Math.Round(p.x, 2), Math.Round(p.y, 2))));
                            }
                            if (p != structEdges[j].b)
                            {
                                //structEdges.Add(new GeometryMath.Segment(p, structEdges[j].b));
                                structEdges.Add(
                                    new GeometryMath.Segment(
                                        new GeometryMath.Point(
                                            Math.Round(p.x, 2), Math.Round(p.y, 2)),
                                        new GeometryMath.Point(
                                            Math.Round(structEdges[j].b.x, 2), Math.Round(structEdges[j].b.y, 2))));
                            }
                            // Удалим старые ребра (i и j)
                            int firstDelete = Math.Min(i, j);
                            int secondDelete = Math.Max(i, j) - 1;
                            structEdges.RemoveAt(firstDelete);
                            structEdges.RemoveAt(secondDelete);
                            // начать с 0;
                            i = -1;
                            j = -1;
                            break;
                        }
                    }
                    j = -1;
                }



                // 3) Удалить повторяющиеся ребра. [нет]
                i = -1; j = -1; // чтоб начать с 0;
                while (i < structEdges.Count - 1)
                {
                    i++;
                    while (j < structEdges.Count - 1)
                    {
                        j++;
                        if (i == j) continue;

                        if (GeometryMath.funcSegment.is_identical(structEdges[i], structEdges[j]))
                        {
                            structEdges.RemoveAt(i);
                            // начать с 0;
                            i = -1;
                            j = -1;
                            break;
                        }
                    }
                    j = -1;
                }

                // Округлим координаты до 2х знаков:
                i = 0; j = 0;
                for (i = 0; i < structEdges.Count; i++)
                {
                    structEdges[i].a.x = Math.Round(structEdges[i].a.x, 2);
                    structEdges[i].a.y = Math.Round(structEdges[i].a.y, 2);
                    structEdges[i].b.x = Math.Round(structEdges[i].b.x, 2);
                    structEdges[i].b.y = Math.Round(structEdges[i].b.y, 2);
                }
                

                // 4) Сформировать список уникальных точек.
                List<Triangula.Vertex> vertexList = new List<Triangula.Vertex>();
                i = 0; j = 0;
                for (i = 0; i < structEdges.Count; i++)
                {
                    // Первая вершина.
                    bool found = false;
                    for (j = 0; j < vertexList.Count; j++)
                    {
                        if ((Math.Round(structEdges[i].a.x, 2) == vertexList[j].x) && (Math.Round(structEdges[i].a.y, 2) == vertexList[j].y))
                        {
                            found = true;
                            break;
                        }
                    }

                    if (found == false)
                        vertexList.Add(new Triangula.Vertex(Math.Round(structEdges[i].a.x, 2), Math.Round(structEdges[i].a.y, 2)));


                    // Вторая вершина.
                    found = false;
                    for (j = 0; j < vertexList.Count; j++)
                    {
                        if ((Math.Round(structEdges[i].b.x, 2) == vertexList[j].x) && (Math.Round(structEdges[i].b.y, 2) == vertexList[j].y))
                        {
                            found = true;
                            break;
                        }
                    }

                    if (found == false)
                        vertexList.Add(new Triangula.Vertex(Math.Round(structEdges[i].b.x, 2), Math.Round(structEdges[i].b.y, 2)));
                }

                // 5) Триангуляция.
                Triangula.Triangulator angulator = new Triangula.Triangulator();
                List<Triangula.Triad> tempTriangles = angulator.Triangulation(vertexList, true);

                // 6) по номерам точек сформировать треугольники List<double[][]>
                List<GeometryMath.Triangle> subTriangles = new List<GeometryMath.Triangle>();
                i = 0; 
                for (i = 0; i < tempTriangles.Count; i++)
                {
                    subTriangles.Add(
                        new GeometryMath.Triangle(
                            new GeometryMath.Point(vertexList[tempTriangles[i].a].x, vertexList[tempTriangles[i].a].y), 
                            new GeometryMath.Point(vertexList[tempTriangles[i].b].x, vertexList[tempTriangles[i].b].y), 
                            new GeometryMath.Point(vertexList[tempTriangles[i].c].x, vertexList[tempTriangles[i].c].y)
                            ));
                }

                // 7) Проверить каждую структурную линию на пересечение с каждым ребром каждого треугольника.
                i = -1; j = -1;
                while (i < structEdges.Count - 1)
                {
                    i++;
                    // точки для линии-разделения.
                    GeometryMath.Point a = new GeometryMath.Point(structEdges[i].a.x, structEdges[i].a.y);
                    GeometryMath.Point b = new GeometryMath.Point(structEdges[i].b.x, structEdges[i].b.y);
                    // список точек для перераспределения.
                    List<GeometryMath.Point> points = new List<GeometryMath.Point>();
                    // 
                    
                    for (j = 0; j < subTriangles.Count; j++)
                    {
                        // Если ребро пересекает треугольник
                        if (GeometryMath.funcTriangle.triangle_intersect(subTriangles[j], structEdges[i].a, structEdges[i].b))
                        {
                            points.Add(subTriangles[j].a);
                            points.Add(subTriangles[j].b);
                            points.Add(subTriangles[j].c);
                            //
                            subTriangles.RemoveAt(j);
                            //
                            j--;
                        }
                    }

                    // удалить повторяющиеся точки.
                    for (j = 0; j < points.Count; j++)
                    {
                        for (k = 0; k < points.Count; k++)
                        {
                            if (j == k) continue;
                            //
                            if (points[j] == points[k])
                            {
                                points.RemoveAt(j);
                                j--;
                                break;
                            }
                        }
                    }

                    //TODO ?: Для каждой точки рассчитать ее перпендикуляр к ребру

                    //TODO ?: Упорядочить полученный массив относительно начальной точки
                    // использовать метрику расстояния между точками.

                    //TODO ?: для каждого участка выполнить процедуру ниже:
                    // Разделим на левую/правую часть.
                    if (points.Count > 0)
                    {
                        List<Triangula.Vertex> leftVertex = new List<Triangula.Vertex>();
                        List<Triangula.Vertex> rightVertex = new List<Triangula.Vertex>();
                        //
                        int onLineCount = 0;
                        //
                        for (int n = 0; n < points.Count; n++)
                        {
                            if (GeometryMath.funcLine.point_in_line_ADVANCED(a, b, points[n]) == 0)
                                onLineCount++;
                        }
                        //
                        if (onLineCount == 2)
                        {
                            for (int n = 0; n < points.Count; n++)
                            {
                                int position = GeometryMath.funcLine.point_in_line_ADVANCED(a, b, points[n]);
                                switch (position)
                                {
                                    case -1:
                                        leftVertex.Add(new Triangula.Vertex(points[n].x, points[n].y));
                                        //points.RemoveAt(n);
                                        break;
                                    case 1:
                                        rightVertex.Add(new Triangula.Vertex(points[n].x, points[n].y));
                                        //points.RemoveAt(n);
                                        break;
                                    case 0:
                                        //points.RemoveAt(n);
                                        break;
                                }
                            }
                            //
                            leftVertex.Add(new Triangula.Vertex(a.x, a.y));
                            leftVertex.Add(new Triangula.Vertex(b.x, b.y));
                            //
                            rightVertex.Add(new Triangula.Vertex(a.x, a.y));
                            rightVertex.Add(new Triangula.Vertex(b.x, b.y));
                            //
                            Triangula.Triangulator angulatorLeft = new Triangula.Triangulator();
                            List<Triangula.Triad> tempTrianglesLeft = angulatorLeft.Triangulation(leftVertex, true);
                            //
                            Triangula.Triangulator angulatorRight = new Triangula.Triangulator();
                            List<Triangula.Triad> tempTrianglesRight = angulatorRight.Triangulation(rightVertex, true);
                            //
                            for (int n = 0; n < tempTrianglesLeft.Count; n++)
                            {
                                subTriangles.Add(
                                    new GeometryMath.Triangle(
                                        new GeometryMath.Point(leftVertex[tempTrianglesLeft[n].a].x, leftVertex[tempTrianglesLeft[n].a].y),
                                        new GeometryMath.Point(leftVertex[tempTrianglesLeft[n].b].x, leftVertex[tempTrianglesLeft[n].b].y),
                                        new GeometryMath.Point(leftVertex[tempTrianglesLeft[n].c].x, leftVertex[tempTrianglesLeft[n].c].y)
                                        ));
                            }
                            //
                            for (int n = 0; n < tempTrianglesRight.Count; n++)
                            {
                                subTriangles.Add(
                                    new GeometryMath.Triangle(
                                        new GeometryMath.Point(rightVertex[tempTrianglesRight[n].a].x, rightVertex[tempTrianglesRight[n].a].y),
                                        new GeometryMath.Point(rightVertex[tempTrianglesRight[n].b].x, rightVertex[tempTrianglesRight[n].b].y),
                                        new GeometryMath.Point(rightVertex[tempTrianglesRight[n].c].x, rightVertex[tempTrianglesRight[n].c].y)
                                        ));
                            }

                            // сместить назад;
                            i = -1;
                        }
                        //
                    }
                }

                // Вывод результатов:
                List<double[][]> resultTriangles = new List<double[][]>();

                i = 0; j = 0; k = 0;
                for (i = 0; i < subTriangles.Count; i++)
                {
                    double[][] tempTriangle = new double[3][];
                    for (j = 0; j < tempTriangle.Length; j++)
                        tempTriangle[j] = new double[2];

                    tempTriangle[0][0] = subTriangles[i].a.x; tempTriangle[0][1] = subTriangles[i].a.y;
                    tempTriangle[1][0] = subTriangles[i].b.x; tempTriangle[1][1] = subTriangles[i].b.y;
                    tempTriangle[2][0] = subTriangles[i].c.x; tempTriangle[2][1] = subTriangles[i].c.y;

                    resultTriangles.Add(tempTriangle);
                }

                return resultTriangles;
            }

            /// <summary>
            /// Принадлежность малых треугольников большим (Списками).
            /// </summary>
            /// <param name="smallT">Малые треугольники.</param>
            /// <param name="bigT">Большие треугольники.</param>
            public static List<List<bool>> checkIsTriangleIn(List<double[][]> smallT, List<double[][]> bigT)
            {
                return GeometryMath.funcTriangle.is_triangle_in_triangle_ADVANCED(smallT, bigT);
            }
        }
    }
}